function [ h ] = MyDeconv2DImage( f,g,type, iterations )
%MYDECONV2DIMAGE Summary of this function goes here
%   Detailed explanation goes here

    if(nargin < 2)
        error('Not enough input arrguments');
    end
    if(nargin < 3)
        type = 'FFT';
    end
    if(nargin < 4)
        iterations = 5;
    end

    if( strcmp(type,'RichardLucy'))
        h = zeros(size(f));        
        for i = 1:size(f,3)
            h(:,:,i) = RLDeconv2D(f(:,:,i),g,iterations);
        end
        return;
    end
    
    if(strcmp(type, 'Bayesian'))
        I = size(f,1)-size(g,1)+1;
        J = size(f,2)-size(g,2)+1;
        h = zeros(I,J,size(f,3));
        for i = 1:size(f,3)
            h(:,:,i) = Bayesian(f(:,:,i),g,I,J,iterations, 0);
        end
        return;
    end
    
    %Type => FFT
    h = zeros(size(f));
    G = fft2(g, size(f,1), size(f,2));
    for i = 1:size(f,3)
        Fi = fft2(f(:,:,i));
        h(:,:,i) = ifft2(Fi./G);
    end
end

function W = RLDeconv2D(H,PSF,iter)

    W = 0.5*ones(size(H));
    
    PSF_HAT = PSF(end:-1:1, end:-1:1);
    
    for i = 1:iter
        x = conv2(W, PSF, 'same');
        x = H./x;
        x = conv2(x, PSF_HAT,  'same');
        W = W.*x;
    end

end

function [ W ] = Bayesian( H,S,I,J, iter, saveing)
    
    [K,L] = size(S);
    SSum = zeros((I+K-1),(J+L-1));
    W = ones(I,J);
    W1 = W;
    for cou = 1:iter
        
        for m = 1:(I+K-1)
            for n = 1:(J+L-1)
                a = max(1,m-K+1);
                b = min(m,I);
                c = max(1,n-L+1);
                d = min(n,J);
                
                x2 = 0;
                for p = a:b
                    for q = c:d
                        x2 = x2+(W(p,q)*S(m-p+1,n-q+1));
                    end
                end
                SSum(m,n) = x2;
            end
        end
        
        for i = 1:I
            for j = 1:J
                y = 0;
                e = i+K-1;
                f = j+L-1;
                
                for m = i:e
                    for n = j:f
                        x1 = H(m,n)*S(m-i+1,n-j+1);
                        y = y+ (x1/SSum(m,n));
                    end
                end
                W1(i,j) = W(i,j)*y;
            end
        end
        W = W1;
        
        if(saveing == 1)
            save(sprintf('tmp/W%d.mat', cou),  'W');
        end
    end
end

